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Flight test and modeling techniques were developed to accurately identify global 
nonlinear aerodynamic models onboard an aircraft. The techniques were developed and 
demonstrated during piloted flight testing of an Aermacchi MB-326M Impala jet aircraft. 
Advanced piloting techniques and nonlinear modeling techniques based on fuzzy logic and 
multivariate orthogonal function methods were implemented with efficient onboard 
calculations and flight operations to achieve real-time maneuver monitoring and analysis, 
and near-real-time global nonlinear aerodynamic modeling and prediction validation testing 
in flight. Results demonstrated that global nonlinear aerodynamic models for a large 
portion of the flight envelope were identified rapidly and accurately using piloted flight test 
maneuvers during a single flight, with the final identified and validated models available 
before the aircraft landed. 
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Nomenclature 

body-axis translational accelerometer measurements, g 
wing span, ft 

wing mean aerodynamic chord, ft 

body-axis nondimensional aerodynamic force coefficients 
body-axis nondimensional aerodynamic moment coefficients 
expected value 

inertia tensor elements, slug-tf 

indicated airspeed, kts 
aircraft mass, slug 

body-axis roll, pitch, and yaw rates, rad/s or deg/s 
number of data points 
dynamic pressure, Ibf/tf 
wing reference area, 

X body-axis component of engine thrust, Ibf 
true airspeed, ft/s 
angle of attack, rad or deg 
sideslip angle, rad or deg 

elevator, aileron, and rudder deflections, rad or deg 
Euler roll, pitch, and yaw angles, rad or deg 
covariance matrix 
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superscripts 


T 

= 

transpose 
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estimate 
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time derivative 

-1 


matrix inverse 


= 

mean 
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normalized 

subscripts 

eg 


center of gravity 

crit 

= 

criteria 

0 

= 

reference value or base term 

P 

= 

prior 

train 

= 

training data 

test 

= 

testing data 


I. Introduction 

A ccurate global aerodynamic models for aircraft are critical for flight dynamic analysis, flight simulation, 
and development of flight control systems. Typically, global aerodynamics are characterized using 
aerodynamic data collected from wind tunnel tests or generated by computational methods such as computational 
fluid dynamics (CFD). These techniques are expensive and time-consuming, and have inherent fidelity limitations 
arising from factors such as wind-tunnel model scale and geometry differences relative to the full-scale aircraft, 
wind-tunnel wall and sting interference, wind-tunnel flow angularity, Reynolds number differences, flow modeling 
deficiencies, and grid geometry approximations for both the aircraft and the flow field. Using flight test methods to 
generate a global aerodynamic model directly avoids all of these problems, although the typical problems associated 
with any flight test still remain, e.g., sensor data quality, achieving adequate data information content from the flight 
test maneuvers, flight test risk and expense, and practical constraints. The advantages of using efficient flight 
testing techniques and nonlinear aerodynamic modeling methods to produce accurate global aerodynamic models 
based on flight data alone include large potential savings in time and cost, in addition to better model fidelity 
resulting from using direct measurements on the real full-scale aircraft in flight, as opposed to using subscale wind- 
tunnel models or mathematical representations of the aircraft and the air flow. In addition, the capability to develop 
high fidelity models in near real time in flight enable a range of possibilities including highly reliable systems and 
self-learning vehicles. 

In previous work, automated multi-axis optimized multisine excitations were applied during time-varying flight 
conditions on a subscale aircraft, and orthogonal polynomial spline modeling techniques were used to identify 
global aerodynamic models from flight data alone Advanced multi-axis piloted maneuvers have been developed 
and flown on an Impala jet aircraft over a large portion of the flight envelope to collect data for global modeling 
using a fuzzy logic modeling technique^. In these works, accurate global aerodynamic models were identified from 
flight data alone, using batch calculations on the ground, after the flight. Other work has examined the problem of 
global aerodynamic modeling from flight data using various modeling methods in post-flight calculations, including 
splines and stepwise regression"^'^, data partitioning with simplified local model structures"^’^’^, multivariate 
orthogonal functions^’"^’^'^^, multivariate B-splines^^’^"^, and combining local simplified models^^'^^. 

The present work is a continuation of the research described in Ref [3], which used advanced piloted 
maneuvers and global aerodynamic modeling based on fuzzy logic, and extends that work to include multivariate 
orthogonal function modeling with splines, onboard data processing, real-time onboard flight data evaluation, and 
global aerodynamic modeling executed during flight test operations. The objective of the flight test project was to 
develop and validate onboard real-time methods for identifying global nonlinear aerodynamic models from piloted 
flight test maneuvers using just a few flight test hours. Flight testing covered a wide range of flight conditions and 
maneuvers, including cruise flight, climbs, descents, turns, stalls, and spins, and novel maneuvers with variations in 
load factor, angular rates, amplitudes, and subsonic Mach number. The integrated real-time global aerodynamic 
modeling approach included real-time onboard monitoring and evaluation of flight data information content, 
efficiently and accurately identifying or updating global aerodynamic models for all six rigid-body degrees of 
freedom, making real-time flight test maneuvering decisions, and conducting model prediction validation tests, all 
onboard the aircraft. 
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This onboard global aerodynamic modeling capability has obvious appeal for saving flight test time and money, 
and for saving post-flight data analysis and modeling time and effort. Perhaps the biggest advantage is real-time 
awareness not only of how well the flight testing is being done at any time, but also of what further maneuvers need 
to be executed to fill gaps in the data collection, using quantitative measures of how well the identified model fits 
the measured data and how well the model predicts aircraft responses. The approach compresses the repeated 
sequence of flight testing, analyzing the data, identifying aerodynamic models, and planning further flight test 
experiments into a single flight (or just a few flights), using onboard real-time methods applied over a wide range of 
flight conditions. 

Real-time onboard global aerodynamic modeling could also be used to identify changes in airplane dynamics, 
enable self-learning vehicles and enhanced control capabilities when merged with advanced control system designs, 
provide information for fault detection and flight envelope protection, and improve flight safety. The capability is 
also a key component in the NASA “Learn-to-Fly” concept, where the intent is to autonomously develop vehicle 
characterization and flight control, with minimal human interaction and analysis time. In addition, the approach 
could be applied to rapidly generate or update the global aerodynamics for flight simulation, based on flight test 
data. The current state-of-the-art in flight testing for global aerodynamic modeling involves processes that are time 
and resource intensive - typically requiring repeated, precisely-flown maneuvers, and significant efforts in post- 
flight analysis. The goal of the current research is to develop flight test and modeling techniques with the potential 
for dramatic reductions in the time and money required to build high-fidelity nonlinear aerodynamic models over a 
large flight envelope. The developed methods execute in real time onboard the aircraft, using flight data alone, and 
produce global aerodynamic models with excellent prediction capability. 

The next section describes the test aircraft. Following that is a description of the piloted flight test maneuver 
technique, and some evaluation and discussion of the resulting flight data. Section IV describes the two global 
aerodynamic modeling techniques applied in flight - one based on fuzzy logic^’^^, and the other based on 
multivariate orthogonal modeling functions^’"^’^’^^’^^. Section V describes how flight test data from the piloted flight 
test maneuvers were used to identify global aerodynamic models, including regions of the flight envelope associated 
with cruise flight, climbs, descents, turns, stalls, and spins. Piloted dynamic maneuvers with power level changes 
were used to include jet engine thrust in the models identified from flight data. Flight tests with various aircraft 
configurations (flap setting and landing gear position) were also included. A flight simulation was created by 
combining previously-developed simulation software"^’^^ written in MATLAB® with the global aerodynamic model 
identified from flight data. This allowed further checks on the prediction capability of the identified global 
aerodynamic model and demonstrated rapid nonlinear simulation development and validation. Details of the flight 
test, real-time modeling, and prediction results appear in Section V. 

All of the real-time flight software used in this work was written in MATLAB® by the authors. Some of the 
software for implementing real-time data analysis and multivariate orthogonal function modeling tasks came from 
the software toolbox called System IDentification Programs for AirCraft, or SIDPAC"^. SIDPAC is bundled with 
Ref [4], and is therefore publicly available^^. 


II. Aircraft 

The test aircraft was an Aermacchi MB-326M Impala, shown in Fig. 1, operated at the National Test Pilot 
School in Mojave CA. The Impala is a 2-seat, tandem single-engine turbojet airplane, used as a trainer or light 
attack aircraft. This particular airplane is fitted with a flight test nose boom equipped with angle-of-attack and 
sideslip angle vanes. Airspeed and dynamic pressure are provided through the production instrumentation system. 
The primary flight control system is completely reversible, with control pushrods connecting the stick and pedals 
with the elevator, ailerons, and rudder. All control surfaces are statically and dynamically balanced. The aileron 
aerodynamic balancing is attained by Irving diaphragms and balance tabs on the trailing edge. The elevator is 
aerodynamically balanced by two tabs on the trailing edge. Trim for roll and pitch forces is implemented though a 
coolie-hat switch on the control sticks. Trim tabs are actuated by electromechanical servos. Trim indications are 
available in both cockpits. No springs or bob weights are included in the control system. Secondary flight controls 
are powered by the aircraft hydraulic system. The speed brake is on the bottom surface of the fuselage, with 
deflection automatically limited to approximately half travel when the landing gear is extended, to prevent the speed 
brake from striking the ground. Flaps are selectable in 3 positions: up, take-off (half), and landing (full). Steering is 
accomplished on the ground by differential braking on the main gear and a castoring nosewheel. The airplane is 
powered by a single Rolls Royce Viper turbojet engine. The engine produces approximately 2500 pounds of thrust 
at sea-level static conditions. Each cockpit is equipped with a Martin-Baker MK- AS. 06 A/M ejection seat that 
allows ejection from the airplane at all speeds and flight altitudes down to zero speed and altitude. A cabin 
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pressurization system maintains cabin pressure and provides heat and air conditioning. A single canopy covers both 
cockpits. 



Figure 1. National Test Pilot School (NTPS) Aermacchi MB 326M Impala aircraft. 

Flight controls and flight instruments are repeated in both cockpits. Flight test maneuvers were executed by the 
pilot in the front cockpit for some flights, but on most flights, the flight test maneuvers were executed by the flight 
engineer / research pilot in the back cockpit. 

Reference geometric characteristics of the airplane are listed in Table 1. A 3 -view drawing of the airplane is 
shown in Fig. 2. The drawing does not include the nose boom fitted onto the test airplane. The longitudinal, 
vertical, and lateral dimensions (FS, WL, and BL, respectively) in Fig. 2 are all in inches. Standard sign convention 
is FS positive aft, WL positive up, BL positive out the right wing. 

The maximum gross weight of the airplane is 9600 pounds (basic airplane with full fuselage fuel and full wing 
tip tanks). For this flight testing, the weight of the fueled airplane with crew was approximately 8185 pounds. The 
allowable CG range is 22% - 30% c . Weight and balance limitations cannot be exceeded by normal operating or 
loading conditions. Mass properties were computed based on measured fuel flow, pre-flight weight and balance, 
and inertia measurements done on the ground for the aircraft with and without fuel. Figure 3 is a diagram showing 
typical CG change with fuel burn for the airplane as tested. Typically the test maneuvers were complete at Bingo 
fuel, as shown in Fig. 3 Table 2 lists the inertia estimates available for the test airplane. 

Measurements that were recorded by the onboard data system and telemetered to the ground included: a, [5, 

Gy, p, q, r, 4 ? 4 ? 4 , q , calibrated airspeed, total temperature, pressure altitude, engine RPM, fuel flow rate, and 
throttle position. Air flow angle vanes installed on the nose boom measured angle of attack and sideslip angle, and 
pressure sensors measured static pressure and total pressure. Potentiometers connected to control system push rods 
were used to measure control surface deflections. Strain gages on the pilot stick and rudder pedal linkages were 
used to measure control forces. Table 3 shows the location of the inertial sensors and air flow angle vanes. A real- 
time analysis computer (RAC) was installed in the aft cockpit (shown in Fig. 4) to conduct analyses and display 
information during flight. Sampling rate for the real-time flight data was 25 Hz. 
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Figure 2. Three-view drawing of the Aermacchi MB 326M Impala aircraft (See Table 1 for dimensions). 
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Figure 3. Longitudinal center of gravity position change with fuel burn 
for the Aermacchi MB 326M Impala aircraft. 
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Figure 4. Aft cockpit of the Aermacchi MB 326M Impala aircraft with RAC installed. 


Table 1. Aircraft Geometry 


c , ft 

6.233 

b , ft 

34.659 

5 ,ft2 

204.5 

Leading edge of c , in 

327.72 

Overall height, ft 

12.20 

Distance between main wheels, ft 

7.78 

Distance from nose gear to main gear, ft 

13.64 

^ ill 

34.59 

Leg ^ in 

0 

z^g , in 

0 
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Table 2. Aircraft Inertia Estimates 


Configuration 

Inertia Tensor Elements, slug-ft^ 

^ XX 

^yy 

tzz 

^xz 

Full Fuselage Fuel, Tip Tanks Full, 
2 Crew 

14515 

7880 

21900 

700 

300 lb Fuselage Fuel, Tip Tanks Empty, 
2 Crew 

4785 

7830 

12075 

700 


Table 3. Aircraft Sensor Locations 


Sensor 

Location, in 

FS 

WL 

BL 

AFIRS (rate gyros and accelerometers) 

334.98 

28 

1.25 

a vane 

124.04 

10 

6.25 

p vane 

120.42 

16.25 

0 


III. Piloted Flight Test Maneuvers 

Flight test maneuver development was an important aspect of this research. Flight data must have sufficient 
information content to enable accurate identification of the global aerodynamics, including nonlinear and unsteady 
aerodynamic effects. Several input types were used during this flight research, including traditional doublet 
sequences and multi-step inputs, and a new type of input called a fuzzy input, which is a piloted quasi-random input 
with varying frequency content and amplitudes. Ah of these inputs were applied by the pilot to implement both 
single-axis and multi-axis excitation. Multi-axis fuzzy inputs were superimposed during gradually changing flight 
conditions (e.g., decelerations and accelerations), in steady turns, in spins, and also through stall and post-stall flight. 
Multi-axis fuzzy inputs are piloted approximations to automated optimized multi-axis excitation inputs that have 
been successfully flight tested on subscale aircraft^ In this work, multi-axis fuzzy inputs were applied to a very 
large flight envelope on a full-scale manned aircraft, including flight conditions representing cruise, maneuvering, 
stalls, spins, departure, and post-stall. 

Any flight test maneuver that used piloted fuzzy inputs was called a fuzzy maneuver. Fuzzy maneuvers 
executed at individual flight conditions were used mainly for model prediction testing. Multi-axis fuzzy inputs 
implemented while the flight condition was varied slowly were used to produce informative flight data for global 
aerodynamic modeling. These global aerodynamic modeling maneuvers were given names like “multi-axis fuzzy 
deceleration maneuver”, etc. 

A. Piloting Technique 

The piloting technique used to execute the flight test maneuvers that generated flight data for global 
aerodynamic modeling will now be described. The general idea is to excite the aircraft dynamic response using 
multi-axis fuzzy inputs with frequency content in a frequency band that encompasses the expected modal 
frequencies for the aircraft dynamic response. At the same time, the nominal flight condition is changed slowly, 
e.g., low angle of attack to high angle of attack, through stall and recovery. Throughout the maneuver, piloted 
inputs for longitudinal stick, lateral stick, and rudder pedals are implemented with the aim of moving the controls in 
an uncorrelated fashion, while covering as much of the entire range of control surface deflections as possible and 
also varying the frequency content of the inputs. This approach builds on previous global modeling maneuver 
design research using automated excitations^ and piloted inputs^. The result is a very high-efficiency global 
maneuver that provides multi-axis flight data with high information content and low correlations among the 
explanatory variables over a large range of flight conditions, for accurate and efficient global aerodynamic 
modeling. 

A comparison of a traditional doublet sequence maneuver at a single flight condition and a multi-axis fuzzy 
deceleration maneuver is shown in Fig. 5. The efficiency of the multi-axis fuzzy deceleration maneuver comes from 
all controls being applied simultaneously throughout the maneuver (as compared to sequentially for the doublet 
maneuver), and from the slowly-varying flight condition, starting from a relatively high speed cruise condition up 
through stall and recovery (as compared to a single flight condition for the doublet maneuver). 
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An even more efficient global flight test maneuver was developed during this work, where occasional power 
level step changes were made during multi-axis fuzzy maneuvers in slowly-changing flight conditions, such as a 
deceleration. This provided information for jet engine thrust modeling, in addition to the data information content 
for global aerodynamic modeling. This type of maneuver was named a powered multi-axis fuzzy deceleration 
maneuver. 

To fly a fuzzy maneuver at an individual flight condition, the pilot initiates the maneuver from a stable trim 
point in straight and level flight. After the trim, the pilot begins the fuzzy maneuver. The pilot excites the aircraft 
dynamic response in pitch, yaw, and roll axes simultaneously. Inputs are varied in direction, magnitude, and 
frequency/speed of movement as much as possible, while maintaining the mean airplane states associated with the 
initial trim. Primarily, the pilot constantly moves the control stick and rudder pedals such that different 
combinations of deflections and rates occur throughout the maneuver, while ensuring that the inputs keep the aircraft 
within load limits, and the mean altitude and airplane attitudes stay close to trim. Depending on the intent of the 
modeling, the inputs can be “large” or “small”. A “large” fuzzy maneuver would attempt to get very large control 
deflections and body rates and accelerations to explore nonlinearities that might be present in those conditions. A 
“small” fuzzy input would involve less extreme control amplitudes and resulting aircraft motion, to produce data 
that could be used by conventional dynamic modeling approaches to identify a local linear model. Using the RAC 
display, the pilot can monitor the real-time data cross-plots of angle of attack and sideslip, and cross-plots of control 
positions, and continue the excitations until these plots are adequately filled in, indicating uncorrelated data covering 
the desired ranges of key explanatory variables. To excite the aircraft adequately, while also obtaining a wide 
variation of control deflections, the fuzzy maneuvers can result in an unusual flight experience, since the objective 
of the inputs is to be uncoordinated. Lateral and longitudinal load factors vary within the planned limits for the test, 
while angular rates and accelerations are also continually varying. With this and all fuzzy maneuvers, sometimes 
the motion coupling with inputs results in some unexpected dynamics, so it is important to be prepared and 
adequately restrained, and to have items secured if planning for the large maneuver inputs that are possible in a 
fighter-class airplane. 

For the multi-axis fuzzy deceleration maneuvers, the pilot again starts at a trimmed flight condition in straight 
and level flight - usually at a reasonable cruise speed to enable data collection from low angles of attack all the way 
up through stall. Similar to the fuzzy maneuver at an individual flight condition, the fuzzy maneuver can be 
conducted using either “large” or “small” inputs for the same purposes as discussed above. When the maneuver is 
begun, the throttle is pulled to idle thrust and the fuzzy piloted inputs are initiated. The fuzzy inputs are similar to 
those used for a single flight condition; however, since the airplane is decelerating, the reference pitch attitude is 
gradually increasing. Additionally, at the higher speeds structural load considerations limit the control deflections 
available to the pilot. To prevent exceeding airplane limits, care must be taken to avoid excessively large control 
inputs until the aircraft decelerates. At lower airspeeds, control effectiveness decreases, requiring larger amplitude 
inputs. These considerations result in larger control surface deflections at the lower airspeeds, and also limit data 
acquisition for high angles of attack at high speeds (limiting Mach effect characterization in the aerodynamic 
modeling). During the deceleration, the pilot executes inputs that have quasi-random amplitude, frequency, and 
phasing for the three axes. This is done continuously during the deceleration and can be reasonably physically and 
mentally demanding since the objective is excitation of the airplane dynamics in an uncoordinated / uncorrelated 
fashion throughout the maneuver. The large amplitude motions used during these maneuvers prevent precise 
altitude control, so to approximate standard level deceleration conditions, the pilot can monitor altitude and adjust 
inputs to remain within a band around the target altitude. To change the rate of deceleration for a maneuver, the 
pilot can select different target values for mean pitch-up attitude and maneuver relative to that. Since air data are 
being measured continuously, adherence to a tight altitude band during the deceleration is not necessary. Because 
the test vehicle was known to have excellent post-stall recovery characteristics, often the deceleration maneuver was 
continued to beyond stall angles of attack. Above stall angles of attack, the airplane dynamics could become larger 
and more unpredictable; however, for modeling purposes, excitation inputs still needed to be made to obtain 
informative data for the aerodynamic modeling in that flight regime. 
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Figure 5. Comparison of conventional doublet sequence inputs and fuzzy inputs. 


Powered multi-axis fuzzy deceleration maneuvers were implemented in the same manner as for the multi-axis 
fuzzy deceleration maneuvers. Initial conditions were a straight and level trim, then the fuzzy control inputs were 
initiated. In addition to the fuzzy stick and rudder pedal inputs described above, for the powered multi-axis fuzzy 
deceleration maneuvers, the throttle was moved in steps from idle to other random power settings during the 
maneuver. Since the engine response is slow relative to the other quantities being varied, best results were obtained 
by setting the throttle and letting the engine speed stabilize before moving the throttle again. To aid in executing the 
maneuver, usually this was done either by selecting airspeed intervals for making a change (for instance, change 
throttle every 15 knots of deceleration) or based on time (e.g., every 10 seconds). Of course, when part of the 
maneuver was being conducted with more than idle power, the time for deceleration to stall was increased, and for 
level flight, increases in power could result in accelerations. Therefore, variation of the deceleration rate for the 
maneuver by setting an initial pitch attitude angle reference (as discussed earlier) worked best. During the 
maneuver, the pilot’s task was to continually vary the input amplitudes, speed of control movement, and phasing of 
pitch, roll, and rudder inputs, while decelerating, maintaining a mean reference attitude, and periodically making 
step changes to the throttle position. All of the same considerations regarding input amplitude applied, as in the 
previous discussions. The maneuver can be conducted until reaching a specified speed, or continued into and 
beyond stall. Engine stability considerations might dictate the cessation of the power steps at high angles of attack, 
but for this airplane, the engine performed well at all flight conditions tested. Figure 6 shows data from a powered 
multi-axis fuzzy deceleration maneuver with an initial pitch attitude reference angle of 20°. This maneuver 
produces very informative data over a large range of nominal angle of attack. Because the nominal angle of attack 
changed slowly, this maneuver can be considered a combination of informative multi-axis excitation data for many 
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different nominal angles of attack, executed in a single, efficient, combined maneuver. This approach is practical 
because the energetic multi-axis fuzzy perturbations excite the aircraft dynamics in a very time-efficient manner to 
produce high data information content, so that the aircraft dynamics can be sufficiently excited even when the flight 
condition is changing. 






Time, s 



Figure 6. Impala flight data for a global aerodynamic modeling maneuver. 


Another maneuver invented for this work was named the Spiral Powered All-axes fuZzy (SPAZ) maneuver. 
This maneuver included a range of angle of attack similar to the fuzzy deceleration maneuvers, but also covered a 
large range of Mach number, load factor, and dynamic pressure during a single maneuver. For this maneuver, a 
block of altitude was used. The maneuver starts at a nominal cruise speed and the nose is pitched over to obtain 
negative angles of attack initially. The uncorrelated all-axes fuzzy inputs are conducted throughout the maneuver as 
explained above - including the step thrust variations. The airplane is put into a descending spiral to build up Mach 
number and dynamic pressure. Gradual increase in load factor is implemented as the speed is increased. 
Eventually, with the added load factor, the nose of the aircraft is raised above the horizon and a deceleration 
commences. This continues - with an upward spiral to enable the variation of load factor without going completely 
vertical - until stall or arrival at the condition for ending the data collection maneuver. 

Executed properly, the SPAZ maneuver includes rich data content across most of the envelope of the airplane in 
one maneuver. Figure 7 shows an example of a SPAZ maneuver. Note that all flight data are continuously varying. 
Mach, dynamic pressure, and mean angle of attack have slow variations, the airplane control surface positions and 
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response variables (instantaneous angles of attack and sideslip, and body-axis angular rates) are all moving at 
varying frequencies throughout the range of the vehicle response capabilities, and the engine RPM (related to thrust) 
is moving in a slower manner that is uncorrelated with the other data. 
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Figure 7. Impala flight data for a SPAZ maneuver. 

Other maneuver types included 3 types of loaded fuzzy maneuvers - fuzzy steady turns, fuzzy wind-up turns, 
and fuzzy wind-up/wind-down turns. Additionally, flights in departed conditions such as spins and spin entries were 
conducted with fuzzy inputs. Most of the previous maneuvers were conducted near 1 g load factor. In order to 
explore Mach effects or other speed-related effects at angles of attack away from trim, maneuvers were needed to 
test at higher load factors. To accomplish this, fuzzy maneuvers were developed for both steady turn conditions and 
for wind-up and wind-down turns. 
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The fuzzy maneuver in a steady turn was conducted by establishing the desired condition (e.g., level 3 g turn at 
250 KIAS). This condition was set as a trim condition, then after established, the fuzzy inputs were initiated. The 
control input strategy was the same as for the other fuzzy maneuvers discussed, except that the reference mean 
conditions used by the pilot were now the initial roll and pitch attitudes, altitude, speed, and load factor. A target 
maximum and minimum load factor were specified for the maneuver, as guidance for the pilot. After the trim, the 
pilot conducted the input about the mean trim conditions until sufficient data were gathered. 





Time, s 


Time, s 


Time, s 


Figure 8. Fuzzy wind-up / wind-down turn. 


Wind-up turns were another maneuver type used to obtain data at higher angles of attack and higher speeds to 
assess Mach or other speed-related effects. The fuzzy wind-up turn started with a level-flight trim condition at a 
specified speed and altitude. The wind-up turn could be conducted either with the addition of power changes or 
with a specified power setting. When the maneuver is started, the mean roll angle and pitch angle is varied to 
maintain a constant mean airspeed while the mean load factor is increased. To reduce pilot workload, after the trim, 
the pilot would roll to approximately 60 degrees bank angle and then begin the fuzzy inputs while increasing load 
factor. As the load factor increased, pitch attitude was lowered to use gravity to augment thrust. The maneuver was 
complete when the target maximum load factor was reached, stall or heavy buffet occurred, or the altitude went 
below a previously specified minimum. 
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The fiizzy wind-up/wind-down turn combined elements of both the fuzzy steady turns and the fuzzy wind-up 
turns. An example is shown in Fig. 8. The maneuver was started with a level-flight trim condition, then began in the 
same way as the fuzzy wind-up turn. When the load factor reached a pre-defmed target, the speed was allowed to 
bleed off while continuing fuzzy inputs and attempting to maintain the target load factor as a mean value. This 
increased the angle of attack and eventually reached stall angle of attack. Typically at that point, the airplane was 
held in the stall/post stall regime as speed continued to bleed, and the end of the maneuver was similar to the fuzzy 
deceleration maneuvers. 

For fuzzy spin or spin entry maneuvers, the airplane was trimmed at a fairly low pre-stall speed (150 KIAS, 
typically). Then, the airplane was decelerated until near stall. At that point, large rudder pedal and aft-stick control 
inputs were made to initiate a departure. During the transient motions, the rudder and pitch inputs were backed off 
slightly, and limited amplitude fuzzy control inputs were made to provide control excitation while maintaining the 
mean control deflections to continue in the departure or spin without rapid recovery. This was done for spin entry 
conditions in both directions, and in developed spins. 

After each maneuver, the flight data were processed onboard the aircraft, to correct for sensor positions and to 
compute necessary quantities such as dynamic pressure, true airspeed, Mach number, nondimensional angular rates, 
and nondimensional aerodynamic force and moment coefficients, which were the measured outputs to be modeled. 
After that, fuzzy logic and multivariate orthogonal function modeling techniques were applied to develop global 
nonlinear aerodynamic models for the aircraft. 

B. Evaluating the Data from the Piloted Flight Test Maneuvers 

With practice, piloted multi-axis fuzzy inputs with a correlation coefficient magnitude of 0.3 or less were 
achieved, indicating a low level of correlation among the inputs for the roll, pitch, and yaw axes. The correlation 
coefficient for two input time series vectors and «2 is computed as 


PuiU2 


N 

Z[“l(0-"l][“2(0-"2] 

i=\ 



( 1 ) 


Although the lowest correlations can be achieved using orthogonal optimized multisine inputs implemented by a 
flight computer^’"^, correlation with magnitude 0.3 is still very low, and well below the threshold value for good 
modeling data, which is 0.9, cf Ref [4]. 

Figure 9 shows cross-plots of measured aircraft states and controls using data from the SPAZ maneuver shown 
in Fig. 7. These plots demonstrate that a wide range of the explanatory variables generally used for aerodynamic 
modeling was swept through during this single maneuver. Note also that the cross-plots generally do not show 
diagonal lines or ellipses, which means that the explanatory variable data from this maneuver had very low 
correlations. Low correlations mean that the aerodynamic dependencies on the explanatory variables can be 
identified accurately and without ambiguity. Cross-plots for other aircraft states and controls were similar in that the 
plots indicated low correlations and wide ranges of coverage for the explanatory variables. 

Multi-axis fuzzy inputs were also flown during accelerations, steady turns, upright spin entries to the left and 
right, climbs, descents, wind-up turns, and other maneuvers. Table 4 contains a summary of what was flown during 
the flight test campaign. 
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Figure 9. Impala flight data cross-plots for a global aerodynamic modeling maneuver. 


IV. Global Aerodynamic Modeling 

The objective of global aerodynamic modeling is to identify a model for each nondimensional aerodynamic force 
and moment coefficient as a function of explanatory variables that can be measured, such as angle of attack, pitch 
rate, and control surface deflections, over a large range of the explanatory variables. There are two main difficulties 
in doing this: 1) experiment design to collect dynamic data over a large range of the explanatory variables, which 
was addressed in Section III, and 2) identifying an accurate global model, which is the subject of this section. 

One approach to global aerodynamic modeling is to identify a series of linear or simplified local models which 
individually are valid for relatively small ranges of the explanatory variables, then join the local models together in 
some way to implement a global aerodynamic model. The fuzzy logic global nonlinear modeling used in this work 
employs a variation of this type of approach, using fuzzy membership functions to combine linear models from 
multiple cells, thereby producing a global nonlinear aerodynamic model^’^^. In another approach, linear and 
nonlinear multivariate modeling functions are orthogonalized to isolate and quantify the modeling capability of each 
individual candidate model term, then a statistical metric computed from the measured data is used to select which 
of the candidate terms should be included in the final model. The multivariate orthogonal function modeling used in 
this work employs this approach to identify compact global nonlinear models valid over large ranges of the 
explanatory variables^’"^’^’^^’^^. For very large ranges of the explanatory variables, or severe local nonlinearities, this 
approach sometimes compromises local model fit in order to achieve a better global model fit. This is the result of 
using global orthogonal polynomial functions for the modeling, which do not have the ability to change locally 
without modifying the entire model. This problem can be remedied by also including spline functions of the 
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explanatory variables in the modeling process, which provides the required local nonlinear modeling capability, 
while retaining clear physical interpretation of the model and the automated selection of optimal modeling 
complexity necessary to accurately characterize the functional dependencies. 

The next subsections describe the general nonlinear aerodynamic modeling problem using nondimensional 
coefficients, model quality metrics used to evaluate the success of the modeling, and the two methods used for 
onboard global aerodynamic modeling in real time. 


A. Nondimensional Aerodynamic Modeling 

For global aerodynamic modeling from flight data, nondimensional aerodynamic force and moment coefficients 
were used as the response variables (also called the dependent variables) for the modeling problem. A separate 
modeling problem was solved for each force or moment coefficient, corresponding to minimizing the equation error 
in each individual equation of motion for the six rigid-body degrees of freedom of the aircraft. Values for the 
nondimensional aerodynamic force and moment coefficients cannot be measured directly in flight, but instead must 
be computed from measured and known quantities using the following equations'^ 


Cx = -C^ - 
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These expressions retain the full rigid-body nonlinear dynamics in the aircraft equations of motion. 

The result is N values of the nondimensional force and moment coefficients, where N is the number of data 
points for a particular maneuver or series of maneuvers. Such data are often called measured force and moment 
coefficient data, even though the data are not measured directly, but rather computed from other measurements and 
known quantities. Data for explanatory variables such as angle of attack, sideslip angle, body-axis angular rates, and 
control surface deflections, are measured directly. 

The desired form of the global aerodynamic model is a mathematical model structure with estimated model 
parameter values and associated uncertainties, relating the nondimensional aerodynamic force and moment 
coefficients to aircraft states and controls that can be measured. All of the global modeling was based on 
equation-error least-squares in the time domain. In this formulation, the dependent variable z, which is one of the 
nondimensional force or moment coefficients, is modeled using an expansion of generally nonlinear model terms 
composed from the explanatory variables, which are nondimensional aircraft states and control surface deflections. 
This leads to an over-determined set of equations for the unknown model parameters, which can be solved using a 
standard least-squares method. 

For example, the equation-error least-squares problem for the pitching moment coefficient could be 
formulated using the model structure 


( 7 ™ — cc + C^ 




+ 5^ + + C^ a5^ + C 






( 7 ) 
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using 


Z — ^0 + £ 


( 8 ) 


where 


: = [C^ (l) (2) ... X 1 vector of values computed from Eq. (4) 


0 = 


x = 


c c c c c 


a adp 1 

2F 


= 5x1 vector of unknown parameters 
= Nx5 matrix of modeling function or regressor vectors 


^ = [^(l) ^(2) ... ^(A^)]^ = A^xl vector of equation errors 


The matrix X is assembled using measured data, with each column representing a modeling function, also 
called a regressor. The best estimator of ^ in a least-squares sense comes from minimizing the sum of squared 
differences between the dependent variable measurements z and the model. 
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The least-squares solution for the unknown parameter vector 6 is' 


0 = 




The estimated parameter covariance matrix is computed from"^ 


X^z 
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The model output is 


y = X0 

and the model fit error variance estimate is 

^2 {z-yf jz-y) 
(N-n) 


( 12 ) 


(13) 


where the number of unknown parameters n = 5 for this example. The standard errors of the estimated parameters 
are given by the square root of the diagonal elements of the covariance matrix. 








(14) 


This approach allows for generally nonlinear modeling functions (such as a 5^ in Eq. (7)), but the solution in 
Eqs. (10)-(13) involves non-iterative linear algebra, which can be done very fast. For global aerodynamic modeling, 
a main issue is identifying which modeling functions should be used for the aerodynamic coefficients, or 
equivalently, determining which functions should appear on the right side of model equations like the example 
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Eq. (7). This important issue is addressed in subsequent sections describing fuzzy logic modeling and multivariate 
orthogonal function modeling. 

B. Modeling Metrics 

Two modeling metrics used to quantify the success of the global aerodynamic modeling. These metrics were 
applied for both fuzzy logic modeling and multivariate orthogonal function modeling. The modeling metrics were: 


2 

7. Coefficient of determination R 


2 ^ M 
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{z-yf (z-y) 
z^z-Nz^ 


(15) 


2 

The R metric quantifies the amount of variation of the dependent variable about its mean value that is explained 

2 2 2 

by the model, so that 0 < 7? < 1 . Often, R is given as a percentage. R is a model fit quality measure. 


2. Predicted Squared Error PSE 
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PSE^ 


[z-xf (z-X0^ 


PSE = 


N 


+ 


N 


2 n 
N TV 


(16) 


(17) 


The PSE metric quantifies expected squared prediction error for an identified model applied to data not used in the 
model identification process. The constant is the upper-bound estimate of the mean squared error between 

future data and the model, i.e., the upper-bound mean squared error for prediction cases. A simple estimate of 

that is independent of the model structure can be obtained by considering to be the residual variance estimate 
for a constant model equal to the mean of the measured dependent variable values. 


Using this simple upper bound estimate for is conservative, and consequently, the PSE metric tends to 

overestimate actual squared prediction errors for new data, cf Ref [21]. Therefore, the PSE metric conservatively 
estimates the squared error for prediction cases. The PSE metric is a combination of a model fit quality measure 
(the first term on the right side of Eq. (16), proportional to the value of the least-squares cost function) and a model 
complexity penalty (the second term on the right side of Eq. (16), proportional to the number of terms in the 
model, n). 

C. Fuzzy Logic Modeling 

The development of the fuzzy logic algorithm was based on previous work ’ ’ , with improvements in the 
underlying membership functions and in several other aspects. The “fuzzification” of the data and the approach 
taken is similar to that proposed originally by Takagi and Sugeno^^. There are two basic approaches to dealing with 
a modeling problem using fuzzy logic. The first uses fuzzy sets that are comprised of membership functions, 
weighting factors, and if-then rules. This process was developed to be similar to the human decision making 
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process. The if-then rules are called linguistic rules and each is single dimensioned (e.g., if the value is high, the 
output is low). Results from this form of fuzzy logic tend to be piece-wise continuous and are used extensively in 
controls applications. The other type of fuzzy logic modeling is based on multi-dimensional fuzzy reasoning. This 
second approach consists of fuzzy cells comprised of internal functions and membership functions. Multiple fuzzy 
cells provide the “fuzziness” of the model, and predicted outputs are smooth. The process involves developing a set 
of fuzzy cells, each cell comprised of linear functions of all of the explanatory variables. A conventional approach 
to dealing with nonlinear modeling is to divide the problem up into a series of piecewise linear solutions. The fuzzy 
approach taken here divides all explanatory variables into ranges, then assigns membership functions to each 
variable. The membership functions for each explanatory variable are combined in all combinations with all the 
other explanatory membership functions. 

7. Normalization of Parameters 

All explanatory variables and output variables used in the fuzzy logic modeling are normalized to values 
between 0 and 1 before starting the modeling process. This is done to prevent numerical biases that could occur 
because of parameters being of different magnitudes, and more importantly, to facilitate the creation and 
interpretation of membership functions for each of the explanatory variables. The constants used for normalization 
need to be retained for transferring the results back from normalized space to real space. The normalization 
procedure is straightforward: 


^hii 


(19) 


where 


^bias 


^min 

^max ~ ^min 


( 20 ) 


The maximum and minimum values can either be computed based on actual values for a given data set, or set a 
priori to cover the range of interest of the explanatory variables. Likewise, the output variable (usually a force or 
moment coefficient) is normalized. For the remainder of the fuzzy model discussion, x will denote the normalized 
variables. 

2. Membership Functions 

A membership function (MF) maps a normalized explanatory variable (e.g. a, P, 5^, etc. ) into ranges of the 
span of the variable. Each explanatory variable can be assigned one or more MF. The MF’s, denoted Afx), have 
values on the interval [0, 1] and vary across the range of the explanatory variable. Values of 0 for the MF equate to 
the assumption that the variable has no effect on the model at that value of the explanatory variable. Conversely, a 
value of 1 for the MF means the associated explanatory variable is very important to the model at that corresponding 
value of the explanatory variable. Using the MF’s is a method of assigning areas in which the variable is important 
to the model and where it is not. 

Membership functions can have many shapes. Traditionally MF’s have been overlapping straight lines, 
triangles, quadratic curves, etc. In the current research, triangles, sinusoidal curves, and ramps have been evaluated. 
All forms had strengths and weaknesses for various applications, however the ramp-shaped MF’s have been used for 
the results shown in this paper, so that is the format that will be discussed. The ramp-shaped membership functions 
are attractive because they give weights to the explanatory variables that have simple physical characteristics. The 
ramp membership functions are designed to have areas of constant I’s and areas of constant O’s with straight lines 
(ramps) connecting the areas. This process divides the explanatory variable range into smaller and smaller sections 
to be used in the modeling as more membership functions are added. 


The ramp membership functions are developed as follows: 


1. Initialize the membership function matrix to 0. The size of the membership function array for an 
explanatory variable is A x Nmp where N is the number of points in the x-variable array, and is the number 
of membership functions for the variable. 


^ = 0 


( 21 ) 
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2. If the total number of membership functions, for a variable x is 1, then: 

A = 1 

3. If there is more than one membership function, for 1 < n < and 1 < i < N: 

a. Set the constant “1” region of the membership function: 

for 

n — 1 n 

then 

Aun = 1 

b. For positive slope ramp region of the membership function: 

for 

72 — 2 72 — 1 

< Xj < 

^mf ^mf 

then 

■)^m/ 

c. For negative slope ramp region of the membership function: 
for 



( 22 ) 


(23) 

(24) 


(25) 


(26) 


then 


N 


< Xi < 


72+1 

N 



(27) 


(28) 


The first membership function for any of the explanatory variables is a constant 1 across the space (Fig. 10). 
This results in a linear fuzzy model, so a linear model can be made by setting all membership function orders to one 
for each variable. If the characteristics being modeled are nonlinear, the parameter space for the explanatory 
variables is divided into progressively smaller segments with higher orders of membership functions. Examples of 
the membership function shapes for a variable with 1, 2, 3, and 6 membership functions are shown in Figs. 10 and 
13. As the number of membership functions increase, the variables are split up into smaller ranges. The number of 
membership functions required for each variable is computed without need for user input during fuzzy model 
structure determination, and will be discussed later in this section. 
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Figure 10. One membership function. 


Figure 12. Three membership functions. 



X 



Figure 11. Two membership functions. 


Figure 13. Six membership functions. 


3. Internal Functions 

Each fuzzy cell is made up of a linear internal function with members from all of the explanatory variables. For 
the current formulation, the format of the internal function is given by: 

= vi + p[xi + + - + plxk (29) 

where the /?’s are constant parameters, y denotes the cell number, and k is the number of explanatory variables being 
used in the model. 

4. Fuzzy Cells 

Numerous fuzzy cells create the fuzziness of the model. Each fuzzy cell is formed by taking one membership 
function from each variable. One fuzzy cell would be: 


fj =pi~\- viA{(x^) + viA{{x2) + ••• + piA^k(Xk) (30) 

where is the membership function from the cell for the explanatory variable. Alternately, the values of the 
membership functions applied to the explanatory variable can be written -, where this would be the product of the 
membership function corresponding to the value with the explanatory variable. The fuzzy model consists 
of the combination of all of the fuzzy cells. The fuzzy cells are created by using all combinations of membership 
functions for each of the explanatory variables - so the final number of fuzzy cells in a model (Nceu) is the product 
of the number of membership functions for each variable. 
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5. Output of Fuzzy Model 

The output of the fuzzy model is defined as the weighted average of the output of the fuzzy cells acting on the 
internal functions. The weight of each cell is assumed to be the product of all membership functions in a given cell. 


yt 




W=fKi*A 


2,i 



(31) 


6. Solving for Fuzzy Model Coefficients 

Using the same equation as above for the output of the fuzzy model and denoting the product of the 
membership functions as 


prod(Aj) = Al^ * A{ i * ... * A{^ (32) 

then the output of the fuzzy model can be re-written as a long linear equation 

1 

* [prod(A^){pl + plx^ i + plx 2 ,i + ••• + +PkXfc,j) (33) 

+ prod(A^)(p^ + + p|x 2 ,i + - + +PfcXfc,i) + - 

+ prod (v4 j + pf‘^‘'“x 2 ,i + — h +Pfc j)] 

The equation for can be multiplied out as a long equation with the many linear parameters (p’s) that can then be 
solved as a linear system using standard least squares parameter estimation techniques. Each cell contribution has 
one parameter to be estimated for each explanatory variable. There are 

Np = N,,ii *k + l (34) 

parameters in total (lumping all the Po terms together) to be identified for the fuzzy model. 

7. Updating the Model Using Recursive Estimation 

As more data becomes available during a test, or if the model fidelity is insufficient, more data can be used to 
update the model. Because this is an equation-error method, additional data could be appended to the existing data 
and the parameter estimation re-done with the new, larger data set. Due to the increasing time required to compute 
the models for progressively longer data arrays, and the large memory that would be required, an alternate approach 
was used. For a fixed model structure, i.e., the selected explanatory variables used for the model and the number of 
membership functions for each of the explanatory variables remain constant, the estimates of the parameters from 
Eq. (33) can be computed using a recursive least-squares approach. 

In the least squares approach, a model for z can be found using the method shown in Eq. (8)-(10). For the 
recursive model updating approach, as developed in Ref [4], updated parameter estimates are arrived at using 

^n+l ~ "f ^n+l[^n+l “ ^n+l^n\ (35) 

where 

^n+l ~ "f -^n+l^n^n+l) (36) 

After new data is added to the model and the parameter estimates are updated, the parameter covariance matrix 
needs to be updated. 


^n+l — 


■ ^n^n+l^n+l^n /(! + 

■^n+l^n-^n+l) 


( 37 ) 
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For the initial parameter estimate, the diagonal elements of the Iq matrix are set to a very large value (with off- 
diagonal elements set to 0), and the initial parameter estimates in Oq are set to 0. With those settings, the final result 
from Eq. (35) is the same as Eq. (10). After each addition of data, the parameter vector and covariance matrix is 
reset for the next model update by 


~ ^n+1 


2:. =2:- 


n+l 


(38) 


8. Determining the Fuzzy Model Structure 

Determining the fuzzy model structure is an iterative process. First, a determination of the approximate 
bandwidth resolution of the model is made. The initial bandwidth assumption is used to select appropriate values 
for thinning of the data, and to compute maximum possible modeling metrics for the measured data used for 
modeling. The data is thinned (downsampled) to reduce computation time. For the data presented here, a cutoff 
frequency of = 2.5 Hz was selected, meaning that only responses at 2.5 Hz or lower were considered. So, the 
thinning factor, TF, is computed 



where / is the sample frequency for the flight data. TF is rounded down to the integer value, then the data is then 
thinned by keeping one point for every TF points in the original array. So, if TF = 1, no thinning occurs. If TF = 2, 
the only half of the points (every other point) are kept. 

The thinned data arrays are then split into training and testing groups. For the splitting, data from the time array 
are alternated at intervals of the 3 times the period of interest, so every 3 /(j0c seconds, the data are alternately selected 
for either the testing or training data groups. 

The next step is to determine the appropriate model structure, which for the fuzzy logic models consists of 
determining the number of membership functions each explanatory variable needs to have. This in turn determines 
the number of fuzzy cells the overall model will have. Following the approach developed in Ref [21], analysis is 
started assuming that all explanatory variables are represented with no membership functions (no explanatory 
variables in the model). Search cycles are conducted by adding one membership function at a time to each of the 
explanatory variables in turn. The number of search cycles is tracked by a variable Ns that is updated at the 
beginning of each search cycle. The coefficients of the internal functions are computed using the training data, and 
the model developed with the training data is then applied to the testing data to obtain prediction of the testing data 
y. The coefficient of determination (F^) for both the testing data fit and the training data model prediction are 
calculated with Eq. (15). Additionally, the predicted squared error {PSE) from Eqs. (16)-(18) is computed after each 
membership function update. 

Each search cycle develops k models to compare (one model for each explanatory variable membership 
function addition). When all explanatory variables have had an additional membership function added and results 
computed, the best model of the search cycle is then determined by picking the combination of membership 
functions that resulted in the best values of a model criterion. Several model criteria were evaluated, and the criteria 
used for these flight tests, F^crit? is the average of the F^ values of the model fit of the training data and the F^ values 
of the prediction for the testing data using the model developed with the training data. The corresponding set of 
membership functions for each variable resulting in the best model based on evaluation of F^crit is then the starting 
point for the next search cycle. After each search cycle, the best model becomes the new basis model (number of 
membership functions for each explanatory variable) and the process is repeated. As more membership functions 
are added, the model fit for the training data generally will continue to improve. When the model becomes over- 
parameterized, the model prediction of the testing data will start to degrade, and at that point the model is complete. 
Additionally, sometimes the addition of membership functions to a variable will slightly improve the model 
prediction, but the improvement is less than the noise in the data or does not improve the overall predictive 
capability of the model. At that point the model would be considered complete. There are 4 criteria, stated below, 
any of which can signal the end of the model determination process. Criteria 1 and 2 were introduced in Ref [23], 
criterion 3 is new to this work. 

Criterion 1: Training correlation coefficient is above the minimum required, and the model criteria correlation 
coefficient is decreasing 
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train (-^5) — min req 

R^cniiNs) > R^cniiK +1) and 

R^UK+^^)>R\itiK+2) 


( 40 ) 


Criterion 2: Training correlation coefficient is above the minimum required, and the model criteria correlation 
coefficient is fluctuating, but generally trending towards decreasing values 

R train(-^5) — R min req and 
R\it(Ns-l)>R\,tmand 

R^cniiNs) < R^cniiK + 1 ) and (41) 

R\nt(N,-l)>R\,i,(N,+l)and 

R^UK-i)>R\,M+2) 

Criterion 3: Training correlation coefficient is above the minimum required, and model criteria correlation 

coefficient is still improving, but predicted squared error is increasing with added cells. For the current modeling, 
trade = 1.10, which means in order to continue to add cells to the fuzzy model, the correlation coefficient 
improvement needs to be at least 10% if the predicted squared error increases by greater than 10%. 

^ trainC^y) — ^ min req ^nd 

^\ri,(V+l)>^WV)and 
R^critiK + 1 ) < trade* R^„itiNs) and 
PSE^UK +1) > PSE„M)/trade and 
R^cniiK +2) < trade* i?^crit(V) and 
f’5£train(V +2) > P5'£'„ai„(V)/trade 


Criterion 4: The model criteria correlation coefficient does not improve after several search updates, or maximum 
number of search updates has been reached. 

After one of the ending criteria is satisfied, the model structure is complete and a final parameter estimation for 
the parameters in all the internal functions of the fuzzy cells is done using all of the data. 

A schematic of the fuzzy logic algorithm and procedures used to analyze flight data is shown in Fig. 14. The 
process steps are outlined and summarized below: 

1 . Data preprocessing: Obtain flight data and make the required corrections to be able to compute overall forces 
and moments acting on the aircraft and the states associated with the flight conditions. This includes such things 
as airspeed and altitude corrections, conversion to pressures, computation of ambient air temperature from 
measured total temperature, sensor position corrections for data from the angle of attack vane, sideslip angle 
vane, and accelerometers, etc. 

2. Selection of Explanatory Variables: The flight variables that the algorithm will use to develop a mathematical 
model of the aerodynamics are selected. Although no functional form of the model is needed (such as 
interactions between variables, etc.), the variables that will be used in a resulting model need to be specified. 
These are selected, along with the computed forces and moments to be used for the identification. The data are 
partitioned into training and testing databases. 

3. Calculation of Membership Functions: Membership functions are used to divide explanatory variables into 
various ranges. The interpretation of membership function values is as follows: where the membership function 
value is 1, the associated explanatory variable is very important in a particular local model over that range of the 
explanatory variable. For a membership function value of 0, the explanatory variable is not important for a 
particular local model over that range of the explanatory variable. When only one membership function is 
selected, it has a constant value of one over the entire range of the explanatory variable. With only one 
membership function for a selected explanatory variable, the fuzzy model is identical to a standard linear model 
relative to that explanatory variable. 

4. Calculation of Internal Function Coefficients for all Fuzzy Cells: The fuzzification of the data is accomplished 
by the use of many fuzzy cells, all comprised of linear models using the explanatory variables through their 
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associated membership functions. Each fuzzy cell internal function includes a combination of one membership 
function for each explanatory variable. Currently, an exhaustive set of all combinations of membership 
functions for each explanatory variable is used to populate the fuzzy cells. Therefore, the models can become 
large. With the sets of internal hmctions, using either Newton gradient descent methods or linear least squares 
procedures to minimize the sum of squared errors, all the coefficients of the internal functions are computed 
using the set of data partitioned for training. 

5. Calculate Model Output: With the fuzzy cells defined, the estimates of the model are computed using a 
weighted output of the fuzzy cells. 

6. Check for Adequacy of Model: The multiple correlation coefficient for the fit of the model to the training data 
is evaluated. If it is above a threshold (indicating a good fit), then that fuzzy model is applied to the set of data 
for testing. If the fit of the training data is not good enough, or if the testing data results continue to show a 
better fit of the data as more membership functions are added, more membership functions are added to the 
model, and the training and testing processes are repeated. Additional checks are conducted to evaluate the 
impact of adding more membership functions relative to overall PSE for the model. If the added membership 
functions do not improve the PSE of the model (i.e. model improvement with more terms is within the noise of 
the available data), and do not improve the model prediction more than a threshold value, then the algorithm 
stops adding membership functions and the model is considered complete. 

7. Compute Final Fuzzy Model: After reaching an acceptable fit to the training data, and when the testing data 
value starts to decrease or when adding more membership functions result in increasing predicted squared error 
for the model predictions, the fit is considered complete and all of the data (training and testing) is used to 
compute the final parameters of all the internal functions in the fuzzy cells. 
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Figure 14. Flow chart for fuzzy logic modeling. 


To make the fuzzy modeling fast enough for real-time flight computations, the model development process 
described above was conducted prior to flight, using data from previous flights to determine the number of 
membership functions for each explanatory variable that were likely to be required to capture the global nonlinear 
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characteristics of the aircraft aerodynamics. This determined the makeup and number of cells required for a model 
that was sufficient to capture the characteristics seen in representative data. Then, for inflight modeling, the pre- 
determined membership function structure was used for the model development. Due to the nature of the fuzzy 
logic algorithm, this might not result in the most efficient model if extra cells were retained in the model, but if the 
cells were not needed for the modeling, the coefficients of those cells would be expected to be very small and would 
not significantly adversely affect the results. By pre-selecting the number of membership functions for each 
explanatory variable, the iterative loop of the process is not needed, and instead a constant overall model structure is 
used as more information becomes available for modeling. This permits additional data to be used to further 
develop and refine the model without having to retain past data. Of course, for initial models of limited flight 
envelope coverage, this approach can lead to over-parameterization issues, degrading predictions. For real-time 
execution of the complete model development process, significant increases of computing power would need to be 
available in the airplane beyond what is currently available with the RAC. 

D. Multivariate Orthogonal Function Modeling with Splines 

The multivariate orthogonal function modeling was based on previous work^’"^’^’^^’^^, with some modifications 
and extensions to implement real-time onboard operation. The technique begins by generating candidate 
multivariate polynomial functions of the selected explanatory variables, up to a selected maximum model 
complexity. Spline functions of the explanatory variables are also included to enable accurate local nonlinear 
modeling, as will be discussed later. Although any function of the explanatory variables could be used, multivariate 
polynomials and spline functions are preferred because of their similarity to a truncated Taylor series and their easy 
physical interpretation. These ordinary functions are then orthogonalized using a Gram-Schmidt numerical 
procedure, so that each of the resulting orthogonal functions retains only the explanatory capability that is unique to 
that modeling function. With this data transformation, it is a straightforward process to select which of the 
orthogonal modeling functions are most effective in modeling the measured data for the dependent variable, and 
how many of these orthogonal functions should be included to identify a model that exhibits both a good fit to the 
measured data and good prediction capability for future data. The final steps are an error-free transformation from 
the selected orthogonal modeling functions back to physically-meaningful polynomial and spline functions of the 
explanatory variables, and calculation of the uncertainties for the associated model parameter estimates. 

7. Generating Orthogonal Modeling Functions 

Multivariate orthogonal functions can be generated from ordinary multivariate functions in the explanatory 
variables using a Gram-Schmidt orthogonalization procedure. This approach is described in Refs. [2], [4], [9], [10], 
and [12], which are the basis for the material presented here. 

The process begins by choosing one of the ordinary multivariate functions as the first orthogonal function. 
Typically, a vector of ones (associated with the bias term in the model) is chosen as the first orthogonal function, 

Pi=l (43) 

In general, any function of the explanatory variables can be chosen as the first orthogonal function, without any 
change in the procedure. To generate the next orthogonal function, an ordinary multivariate function is made 
orthogonal to the preceding orthogonal function(s). Define the orthogonal function pj as 

7-1 

Pj = - Z ykj Pk J = 2, 3,..., itf (44) 

k=l 

where is the ordinary multivariate function vector, and n^ is the total number of ordinary multivariate 
functions to be orthogonalized. Each could be any function of the explanatory variables, such as a multivariate 

polynomial function or a spline function. The for k = \,2, ...,j are scalars determined by multiplying both 

sides of Eq. (44) by , then invoking the mutual orthogonality of the pj^, k = 1,2, j\ and solving for 
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(45) 


Ykj 


pUj 

Pkpk 


The same process can be implemented in sequence for each ordinary multivariate function j = 2,3, . 

The total number of ordinary multivariate functions used as raw material for generating the multivariate orthogonal 
functions, including the bias term, is . Equations (43) - (45) show that each orthogonal function can be expressed 

exactly in terms of a linear expansion of the original multivariate functions. The orthogonal functions are generated 
sequentially by orthogonalizing the original multivariate functions with respect to the orthogonal functions already 
computed, so that each orthogonal function can be considered an orthogonalized version of an original multivariate 
function. 

If the pj vectors and the vectors are arranged as columns of matrices P and X , respectively, and the 
are elements in the row and column of an upper triangular matrix G with ones on the diagonal. 


Then 


G = 


1 Yn 

0 1 
0 0 


Yu "• Yin, 
Y23 "• Yin, 
^ " Yin, 


0 0 0 ••• 1 


X = PG 


which leads to 


P = XG~^ 


(46) 


(47) 


(48) 


The columns of G~^ contain the coefficients for expansion of each column of P (i.e., each multivariate 
orthogonal function) in terms of an exact linear expansion in the original multivariate functions in the columns of 
X . Equation (48) can be used to express each multivariate orthogonal function in terms of the original multivariate 
functions. The manner in which the orthogonal functions are generated allows them to be decomposed without 
ambiguity into an expansion of the original multivariate functions, which have physical meaning. 

2. Nonlinear Identification using Orthogonal Modeling Functions 
The form of a multivariate orthogonal function model is 

z = a^p^+a2P2 + ... + %Pn+s 

where z is an A^-dimensional vector of dependent variable data (e.g., nondimensional force or moment coefficient), 
Z = , modeled in terms of a linear combination of n mutually orthogonal modeling functions 

Pj, j = l,2,...,n . Each pj is an A^-dimensional vector that in general depends on the explanatory variables. The 
aj, j = 1,2,..., n are constant model parameters to be determined, and ^denotes the modeling error vector. 

Equation (49) represents a mathematical model for the dependent variable z in terms of orthogonal functions 
generated from the explanatory variable data. The important question of determining which modeling functions 
should be included in Eq. (49), which implicitly determines n , will be addressed now. 

Assembling the n orthogonal modeling functions from Eq. (49) in the columns of an Nxn matrix P , 

P = [Pl>P2’ ->Pn] ( 50 ) 


26 

American Institute of Aeronautics and Astronautics 



and defining the unknown parameter vector a = ]^ , Eq. (49) can be written as 

z = Pa-\-£ (51) 

Eq. (51) is the same model equation discussed earlier (cf Eq. (8)), except that the modeling functions are now 
orthogonal functions. When the modeling functions are multivariate orthogonal functions generated from the 
explanatory variable data, it is easier to determine an appropriate model structure, because the explanatory capability 
of each modeling function is completely distinct from all the others. This decouples the least squares modeling 
problem, as will be shown now. 

For mutually orthogonal modeling functions, 

p[pj=0 , i^j , = (52) 


T 

and P P is a diagonal matrix with the inner product of the orthogonal functions on the main diagonal. Using 
Eqs. (50) and (52) in the least-squares solution from Eq.(lO), the element of the estimated parameter vector a is 
given by 


^j=[p] ^)l[p] pj] 

The least-squares cost function using orthogonal modeling functions is then"^’^ 


j(«) 




(53) 


(54) 


Equation (54) shows that when the modeling functions are orthogonal, the reduction in the least-squares cost 
function resulting from including the term aj pj in the model depends only on the dependent variable data z and the 

added orthogonal modeling function pj . The least-squares modeling problem is therefore decoupled, which means 

each orthogonal modeling function can be evaluated independently in terms of its ability to reduce the least-squares 
model fit to the data, regardless of which other orthogonal modeling functions are already selected for the model. 
When the modeling functions pj are instead polynomials in the explanatory variables (or any other non-orthogonal 

function set), the least squares problem is coupled, and iterative analysis is required to find the subset of modeling 
functions for an adequate model structure. 

The orthogonal modeling functions to be included in the model are chosen to minimize predicted squared error 
metric PSE defined in Eq. (16) 


J\[ J\[ 


(55) 


The PSE depends on the mean squared fit error, (^z~PdY i^z-Pd^jN = 2 , and a term proportional to 

the number of terms in the model, n . The latter term prevents over-fitting the data with too many model terms, 
which is detrimental to model prediction accuracy'^’^^ While the mean squared fit error must decrease with the 

addition of each orthogonal modeling hmction to the model (by Eq.(54)), the over-fit penalty term cj^ax ^1^ 
increase with each added model term (n increases). Introducing the orthogonal modeling functions into the model in 

order of most effective to least effective in reducing the mean squared fit error (quantified by /?J Pj j for 

the 7 ‘th orthogonal modeling function) means that the PSE metric will always have a single global minimum. 
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Figure 15 depicts this graphically, using actual modeling results from Ref. [4]. The figure shows that after the 
first 6 modeling functions, the added model complexity associated with an additional orthogonal modeling function 
is not justified by the associated reduction in mean squared fit error. This point is marked by a minimum PSE, 
which defines an adequate model structure with good predictive capability. Ref [21] contains further statistical 
arguments and analysis for the PSE metric, including justification for its use in modeling problems. 


X 10® 



Figure 15. Model structure determination using orthogonal functions and PSE. 

Using orthogonal functions to model the dependent variable data makes it possible to evaluate the merit of 
including each modeling function individually, based on the PSE metric. The goal is to select a model structure with 
minimum PSE, and the PSE always has a single global minimum for orthogonal modeling functions. This makes 
the model structure determination a well-defined and straightforward process that can be (and was) automated. 

3. Conversion to Physically-Meaningful Multivariate Function Models 

After the model structure is determined using multivariate orthogonal modeling functions for minimum PSE, 
the identified model output is 


y = Pd (56) 

where the P matrix now includes only the n orthogonal functions selected in the model structure determination, 
n<n,. Each retained orthogonal modeling function can be decomposed without error into an expansion of the 

original multivariate functions in the explanatory variables, using the columns of G~^ in Eq. (48) corresponding to 
the retained orthogonal functions. Common terms are combined using double precision arithmetic to arrive finally 
at a model using only original multivariate functions in the explanatory variables. Terms that contribute less than 
0.1 percent of the final model root-mean- square magnitude are dropped. 

Because each selected orthogonal modeling function can be decomposed without error into a linear combination 
of the original multivariate functions, the variances computed for parameters in the orthogonal function parameter 
vector estimate d (diagonal elements of the 2^ («) matrix) can be combined to compute the variances for the model 

parameters associated with the modeling functions in the final model. Standard errors for the final model 
parameters are then obtained using Eq. (14). 

The final form of the model is a sum of ordinary multivariate functions in the explanatory variables, with 
associated model parameter estimates and standard errors. An example of the final model form obtained is given in 
the next section. 

4. Including Spline Functions in the Ortho gonalization 

For global aerodynamic modeling, functional dependencies of the aerodynamic coefficients on the explanatory 
variables can exhibit significant localized variation. In these cases, a global polynomial model can be inadequate for 
capturing the local variations. To solve this problem, spline functions in the explanatory variables can be introduced 
as additional pseudo-explanatory variables. Splines have local modeling capability that global polynomial modeling 
functions do not have, because splines are polynomials defined only on selected intervals. Low-order polynomial 
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terms defined on limited intervals can approximate nonlinearities quite well. Using polynomial splines also retains 
clear physical interpretation in the final model. 

Spline functions are defined as piecewise polynomials functions of degree m in one or more explanatory 
variables. The term “piecewise” means that the polynomial is different for specific ranges of the explanatory 
variables. Spline function values and derivatives agree at the points where the piecewise polynomials join. These 
points are called knots, and are defined as specific values of each explanatory variable. A polynomial spline (v) 

of degree m with continuous derivatives up to degree m -1 , for a single explanatory variable x e 
expressed as 


where 






/ xm \{x-xX 


X > x^- 
X < x^- 


(57) 


(58) 


and the are constants. The individual (x - x^ terms are the piecewise parts of the spline function. The values 
Xj , X 2 , . . . , Xy^ are knots that satisfy the condition 


Xo <Xi <X 2 <...<x^ <x^^^ (59) 

Three piecewise spline types are sketched in Fig. 16 for the same selected knots in angle of attack. 


(a) 


D liO 20 30 ^0 



0 10 20 30 ^0 

& (deg) 



a id&gj 

Figure 16. Polynomial splines in angle of attack: (a) zero degree, (b) first degree, (c) second degree. 

Note that when the spline knots are the same, a higher-order piecewise spline in a single explanatory variable 
can be computed as a multiplication of lower-order piecewise splines in the same explanatory variable. 
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(60) 


This fact can be used advantageously in generating multivariate orthogonal modeling functions with excellent 
local modeling capability. This is best demonstrated with a simple example. 

Suppose that the nondimensional vertical aerodynamic force coefficient is being modeled with explanatory 
variables angle of attack a and elevator control deflection S^. For a maximum selected model order 2, the ordinary 
polynomial modeling functions that will serve as raw material for the orthogonalization process are: 

l,a,a^ (61) 


which would lead to a multivariate polynomial model of the form 


Cz - Cz„ + Cz„ « + Cz^2 « + Cz„, + Czg + ^Z . 


(62) 


where the values of the model parameters such as and would be estimated from the data using the 

orthogonal function modeling and subsequent decomposition procedures described earlier, and some of the terms 
might not be present, depending on the results of the model structure determination using orthogonal functions. This 
represents a global polynomial model. 

Now introduce a first-order spline term in angle of attack with a single knot located at 10 deg. Using the same 
maximum model order 2, the set of ordinary polynomial modeling functions that will serve as raw material for the 
orthogonalization process expands to 


1 , a, a^, aS^, S^, 


a{a-\0)\,5^{a-\0)\,{a-\0)\ 


(63) 


which would lead to a final multivariate polynomial model of the form 


^Z - + ^Z„“ + ^Z 7 “ + ^Z„* + Q 


-Q,(«-10)UCz ,«(«-10)UCz , <5,(a-10)UCz, («-10) 


(64) 


As before, some of the terms in Eq. (64) might be not be present, depending on the results of the model structure 
determination using orthogonal functions. 

With this expanded set of ordinary polynomial modeling functions, more orthogonal functions will be 
generated, and with the capability to accommodate changes in the linear dependence of on a for a >10 deg , 

as well as additional and aS^ nonlinearities that only take effect when a >10 deg . Note also that a(a-lO)^ 

and (cir-lO)^, for example, are in general not the same function. Because higher-order splines are created 

automatically by the multiplication of first-order splines (cf Eq. (60)), only first-order splines in the explanatory 
variables need to be supplied to the modeling algorithm as pseudo-explanatory variables. 

Extrapolating from this simple example, it is clear that multiple knots in the explanatory variables would provide 
a very effective nonlinear modeling capability in multiple dimensions with local nonlinear modeling capability, 
while retaining physical insight into the functional dependencies. The automated orthogonalization and sorting 
process described earlier identifies which nonlinear terms are necessary to characterize the functional dependencies 
and estimates the associated model parameters. Inputs required from the analyst relate only to the limits of what 
should be considered, such as which explanatory variables to consider, maximum model order to consider, and knot 
locations to consider. However, these can be specified very generously, because the orthogonal function modeling 
algorithm automatically sorts out which of the terms are important, based on the data, and discards the rest. The 
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result is a global parsimonious model with excellent local nonlinear modeling capability and clear physical 
interpretation. 

5. Bayesian Updates 

When the model structure is judged to be adequate, a common task is to update the model parameter estimates in 
that model structure, using data from additional flight test maneuvers. This is useful when individual maneuvers 
have good information content for estimating some parameters but not others, or for incorporating data information 
from more than one maneuver into the model. A Bayesian form of equation-error parameter estimation can be used 
to combine previous modeling results with flight data from an additional maneuver. This approach was used to 
update the global aerodynamic model using multivariate polynomials and splines. The implicit assumptions made in 
doing this, in addition to the usual equation-error assumptions, are: 1) the noise is of the same character for all flight 
data, and 2) the experimental conditions of the flight data are sufficiently similar that the model parameters can be 
expected to be the same. Violating the first assumption leads to inaccuracies in the estimated parameter 
uncertainties. If the second assumption is violated, the parameter estimation results are some weighted average of 
the results that would have been obtained from each subset of flight data analyzed individually. 

To incorporate results from modeling of prior data, the least-squares cost function to be minimized is modified 
as follows"^’^"^ 




(65) 


where Op and Up are from modeling prior data, and the matrix X and vector z are assembled using data from the 

additional maneuver, using the same model structure that was used in the previous modeling. The parameter 
estimation results in this case are obtained from'^’^'^ 




The covariance matrix associated with the parameter estimates is 




-1 


where a is estimated from 


(N-n) 


( 66 ) 


(67) 


( 68 ) 


and n is the number of model terms, equal to the number of columns in the matrix X. 

When the assumptions listed above are acceptable, using prior values for the parameter estimates and covariance 
matrix is equivalent to concatenating the data from more than one maneuver. Using prior values for parameter 
estimates and the covariance matrix is much more convenient and practical than concatenating data, especially for 
onboard calculations involving many maneuvers, as in this work. This approach provides an easy and rigorous way 
to combine the information from additional maneuvers with model results already calculated from prior data. The 
Bayesian update expressions above can be made more general if the noise sequences for the maneuvers to be 
combined have different variances (cf Ref [24]); however, for this work, the expressions given in Eqs. (66)-(68) 
were used. 
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V. Flight Test Results 

Flight testing was conducted in cooperation with the National Test Pilot School in the instrumented Impala 
aircraft, shown in Fig. 1 and described in detail in Section II. 

Flight data information content and flight test maneuver quality could be monitored in real time using the tablet 
computer installed in the aft cockpit. Figs. 17 and 18 are screen shots showing typical displays during flight test 
maneuvers. The displays show real-time strip chart plots and cross-plots of the important explanatory variables. In 
particular, all cross-plots of the control surface deflections were displayed on the screen shown in Fig. 18, to 
evaluate and guide the fuzzy maneuver execution. The real-time data shown in the strip chart plots and cross-plots 
was selectable by the flight engineer in the aft cockpit. Pushbuttons on all screens were used to navigate between 
various screens in the real-time software, and to control the flight data recording. 



Figure 17. Data monitoring screen in the Impala real-time software. 



Figure 18. Maneuver monitoring screen in the Impala real-time software. 
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After the flight test maneuver data was monitored and recorded, the Model pushbutton was selected to switch to 
the modeling screen shown in Fig. 19. This screen was used to evaluate, identify, and/or update the global 
aerodynamic models. When this screen is selected, the real-time software evaluates the current onboard global 
aerodynamic models by applying those models to predict the current flight data, which is the data that has just been 
recorded. Current model performance is evaluated using the and PSE modeling metrics described earlier. If no 
current model exists (i.e., the current data is from the first flight test maneuver), then global models are identified 
based on the current flight data only. Green and red indicator lights on the right side of the screen indicate good 
(green) or inadequate (red) performance of the current model when used to predict the current flight data. Individual 
models or models for all six rigid-body degrees of freedom can be updated using the current flight data by selecting 
pushbuttons on the screen. This implements global model identification for the initial flight test maneuver, and 
sequential model updating for subsequent flight test maneuvers. Fuzzy logic or orthogonal function modeling can 
be selected by a pushbutton toggle on the left side, and the modeling can be re-initialized (i.e., identify new models 
based only on the current flight data) using the red pushbutton on the lower right. Figure 20 shows modeling results 
using multivariate orthogonal function (MOF) modeling, for the same flight data as in Fig. 19 for fuzzy logic (OFL) 
modeling. Transition from Figs. 19 to 20 was achieved by simply selecting the modeling method toggle pushbutton 
on the left side of the screen. 

The modeling metric shown on the right in Figs. 19 and 20 was based on the deterministic portion (i.e., 
excluding random noise) of the measured outputs, computed using a global optimal Fourier smoother"^’^^. If the 
current model characterized 75 percent of all deterministic variation in the measured output (one of the 
nondimensional force or moment coefficients), then the associated model quality indicator was colored green; 
otherwise, the model quality indicator was colored red. The prediction error metric was based on how well the 
current model predicts the current data, compared to the prediction error metric computed from Eq. (16) when the 
model was identified or last updated. If the root-mean- square of the actual prediction error using the current model 
was less than 125 percent of the square root of PSE for the current model, then the associated model quality 
indicator was colored green; otherwise, the model quality indicator was colored red. Model quality indicators were 
shown for the global aerodynamic models identified for each of the six rigid-body degrees of freedom, and for each 
of the two modeling metrics. The green and red lights were used to enable quick evaluation of the modeling results 
in flight. 

Figures 19 and 20 show typical global nonlinear modeling results, based on a multi-axis fuzzy deceleration 
maneuver. Note that all of the model quality indicators are green for both OFL and MOF modeling methods, 
showing that both modeling methods achieved accurate modeling results. 

Throughout the flight testing, both fuzzy logic and multivariate orthogonal function modeling were applied in 
flight. This could be done simply in the screens shown in Figs. 19 and 20, using the modeling method toggle 
pushbutton on the left side. In some individual cases, one or the other method exhibited slightly superior 
performance, but overall, both methods produced similar excellent global nonlinear modeling and prediction results. 
Modeling results were automatically saved when the Save pushbutton was activated, or when the model screen was 
exited, e.g., to go back to one of the screens shown in Figs. 17 and 18 to monitor and record data for the next flight 
test maneuver. 

Jet engine thrust modeling and thrust-induced aerodynamic effects were combined with the global modeling for 
aerodynamic Cx , by modeling the total force along the X body-axis. Ex , and including engine speed in rpm as an 
additional explanatory variable for the global Ex model identification. This generally worked well, as shown in 
Figs. 19 and 20, and avoided the need for a separate model for jet engine thrust. The modeling was done 
dimensionally because the nondimensionalization normally used for aerodynamics (cf Eqs. (2)-(6)) does not apply 
for jet engine thrust. Accelerometer measurement includes effects of both thrust and aerodynamics (cf Eq. (2)), 
so these effects were modeled together in flight. Note that the Ex plots in the modeling screen showed Ex in 
thousands of pounds /lOOO) , for better plot scaling and visual interpretation, cf Figs. 19 and 20. 
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Figure 19. Fuzzy logic (OFL) modeling screen in the Impala real-time software. 



Figure 20. Multivariate orthogonal function (MOF) modeling screen in the Impala real-time software. 
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Figure 21 shows an example of the View screen in the real-time software, which enabled a comparison of the 
identified global nonlinear model with flight data from the current maneuver, using 3-dimensioinal plots with 
selectable plot variables. This was useful for obtaining a global view of the identified model fit to the flight data. 



Figure 21. Identified model 3D View screen in the Impala real-time software. 


Global aerodynamic models could be refined and updated with flight data from more than one maneuver, using 
recursive parameter estimation for fuzzy logic modeling, and the Bayesian approach for multivariate orthogonal 
function modeling, as described earlier. The modeling shown in Figs. 19 and 20 was for the first global modeling 
maneuver on flight 9, a multi-axis fuzzy deceleration maneuver. Following that, a second multi-axis fuzzy 
deceleration maneuver was flown, with prediction results shown in Figs. 22 and 23 for OFL and MOF, respectively, 
showing unsatisfactory predictions. At this point, the Identify All pushbutton was selected, which updated the 
global models, producing the updated modeling results in Figs. 24 and 25 for OFL and MOF, respectively. The 
resulting updated global models were then based on flight data from both maneuvers together, because of the 
sequential analytical methods used to update the global models. Such updates generally improve the quality of the 
identified global models, because of the additional data information content in multiple maneuvers. Subsequent 
predictions for local multi-axis fuzzy maneuvers at 200 KIAS, 150 KIAS, and 100 KIAS are shown in Figs. 26 
through 31 for the OFL and MOF modeling. Except for some slight mismatches in the modeling, very 
successful prediction results were achieved with no further model updates, demonstrating that the global 
aerodynamic modeling from both methods, based on only two global modeling maneuvers, accurately predicted 
local maneuvers over a large portion of the flight envelope. Predictions using OFL modeling could be improved to 
all green metrics if the data from the local fuzzy maneuvers was used to update the Fx model. For MOF, this 
procedure made the model predictions for Fx worse. This is likely the result of the OFL modeling having localized 
cells, and therefore being able to more accurately model localized variations without impacting other parts of the 
model, whereas that was not possible with MOF modeling, which uses mostly global modeling functions and fewer 
model parameters. 

On flight 9, this process was carried out using both fuzzy logic and multivariate orthogonal function modeling 
techniques with similar excellent success for three aircraft configurations: cruise configuration with flaps and 
landing gear up (examples shown in Figs. 19 through 31), take-off configuration with flaps half and landing gear up, 
and approach/landing configuration with flaps full and landing gear down. All of this was done onboard the aircraft 
during a single flight lasting less than 1 hour. A total of 36 global aerodynamic models were successfully identified 
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and validated during this single flight - 6 degrees of freedom, for 3 aircraft configurations, using 2 global modeling 
techniques. All global nonlinear modeling was completed, validated, and saved before the aircraft landed. 



Figure 22. Second global maneuver prediction using OFL. 
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Figure 23. Second global maneuver prediction using MOF. 
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Figure 24. Second global maneuver prediction using updated OFL. 



Figure 25. Second global maneuver prediction using updated MOF. 
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Figure 26. Local fuzzy maneuver prediction at 100 KIAS using OFL. 



Figure 27. Local fuzzy maneuver prediction at 100 KIAS using MOF. 
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Figure 28. Local fuzzy maneuver prediction at 150 KIAS using OFL. 



Figure 29. Local maneuver prediction at 150 KIAS using MOF. 
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Figure 30. Local fuzzy maneuver prediction at 200 KIAS using OFL. 



Figure 31. Local fuzzy maneuver prediction at 200 KIAS using MOF. 
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Table 4 . Impala Flight Data Summary 


Flight 

Primary Maneuvers and Objectives 

1 

Calibrations, system checkout, trim shots 

2 

Multi-axis fuzzy maneuvers at individual flight conditions and during deceleration 

3 

Multi-axis fuzzy, fuzzy loaded deceleration, powered multi-axis fuzzy deceleration maneuvers 

4 

Powered multi-axis fuzzy deceleration maneuvers, departure conditions, fuzzy wind-down turn, 
spin entries 

5 

Powered multi-axis fuzzy deceleration maneuvers, fuzzy wind-up turn, fuzzy wind-down turn 

6 

Powered multi-axis fuzzy deceleration maneuvers, fuzzy wind-up turn, fuzzy wind-up/wind- 
down turn 

7 

Powered multi-axis fuzzy deceleration maneuvers and model validation maneuvers. Clean and 
T/0 configurations 

8 

Powered multi-axis fuzzy deceleration maneuvers, fuzzy wind-up/wind-down turns 

9 

Powered multi-axis fuzzy decelerations and model validation maneuvers. Clean, T/0, and 
Landing configurations 

10 

Powered multi-axis fuzzy deceleration, SPAZ, fuzzy spin entry, and model validation maneuvers 

11 

SPAZ, model check, aerobatic maneuvers, powered multi-axis decelerations and model 
validation maneuvers in landing configuration 

12 

Full-envelope identification with SPAZ maneuvers, powered multi-axis fuzzy deceleration, 
fiizzy wind-up/wind-down turn, fuzzy spin entry and model validation maneuvers. Clean and 
with speed brake extended 


Flight 12 was intended to examine the capability for identifying a global nonlinear aerodynamic model from 
expanded flight envelope data, including large Mach and load factor variations, large ranges of angle of attack and 
sideslip angle, post-stall flight, and spin entries, ah in a single flight. Several maneuvers were flown to gather 
information for model development, then several quick checks across the flight envelope were conducted to verify 
the model with different maneuvers. Individual maneuvers are listed in Table 5. Ah maneuvers included random 
step changes in power level. Flight data from the first maneuver (fuzzy deceleration. Card 3) was used to identify 
the initial global nonlinear aerodynamic model. The next maneuver (SPAZ maneuver. Card 4) was a more dynamic 
maneuver encompassing a larger portion of the flight envelope. The initial model was used to predict the flight data 
from the SPAZ maneuver, and any coefficients with model metrics not meeting the desired values (red indicators) 
were updated using this new data. Five coefficients were updated for both the fuzzy logic (OFL) and multivariate 
orthogonal function (MOF) techniques, using the individual Identify pushbuttons. The resultant model was then 
used to predict flight data from the following maneuver, which was another fuzzy deceleration maneuver. Card 5. 


Table 5. Flight 12 Model Identification Maneuvers 


Card # 

Maneuver Type 

Description 

3 

Fuzzy Deed 

Initial conditions 250 KIAS at 15K Hp. Mean climb angle of 20 degrees during 
deed with fuzzy control inputs 

4 

SPAZ 

Initial conditions 200 KIAS at 20K Hp. Right turning spirals up to 5-g with fuzzy 
control inputs, decreasing speed at conclusion to stall at 50°-60° nose high attitude 

5 

Fuzzy Deed 

Initial conditions 250 KIAS at 15K Hp. Mean climb angle of 30 degrees during 
deed with fuzzy control inputs 

6 

Fuzzy Deed 

Initial conditions 250 KIAS at 15K Hp. Mean climb angle of 60 degrees during 
deed with fuzzy control inputs 

7 

WUT/WDT 

Initial conditions 220 KIAS at 20K Hp. WUT to 4g, then deed in buffet hold to 
stall, with fuzzy control inputs 

8 

Fuzzy Right Spin 
Entry 

Initial conditions 150 KIAS at 18K Hp. Fuzzy inputs to departure with rudder 
bias for left departure. Recover after one turn. 

9 

Fuzzy Left Spin 
Entry 

Initial conditions 150 KIAS at 18K Hp. Fuzzy inputs to departure with rudder 
bias for right departure. Recover after one turn. 

13 

SPAZ 

Initial conditions 200 KIAS at 20K Hp. Left turning spirals up to 5-g with fuzzy 
control inputs, decreasing speed at conclusion to stall at 50°-60° nose high attitude 
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The model prediction of the Card 5 fuzzy deceleration maneuver data is shown in Figs. 32 and 33. The 
predictions are generally good for both techniques as can be seen by the match of the time histories of the flight data 
and the model values. For OFL, all the model fit metrics (R^) were green; however, none of the predicted error 
metrics were green. Therefore, all coefficients were updated with the new Card 5 data included, resulting in a model 
with all modeling metrics green. The MOF method results were similar, with all but one model fit metric green, and 
all but one prediction error metric red. So, all coefficients were updated to include the new Card 5 data, resulting in 
a model with all modeling metrics green. Note that the modeling metrics in Figs. 32 and 33 are lowered because of 
high-frequency buffeting (visible in the flight data for the aerodynamic coefficients), which was not captured by the 
modeling, because it is a high-frequency unsteady aerodynamic effect. This unmodeled aerodynamic effect was 
therefore classified as noise for this modeling, which adversely affected the prediction error metric, and to a lesser 
extent, the model fit metric. This was reflected in the model metric indicators for this prediction case, and others 
like it that involved flight at angles of attack near stall and post-stall, where buffeting occurred. 

The new model was used to predict the data from the next maneuver, another fuzzy deceleration, but at a very 
steep pitch angle (Card 6). The model predictions for the OFL and MOF models are shown in Figs. 34 and 35. The 
models identified with both techniques produced excellent predictions of the global aerodynamics computed from 
the flight data for this maneuver. No changes to either the OFL or MOF models were made at this point since all the 
model metrics were green. Note that Figs. 34 and 35 are model predictions, because the Card 6 data has not been 
used in the modeling. Comparing Figs. 32 and 33 (an earlier prediction case) with Figs. 34 and 35 shows that the 
global modeling improved and the global aerodynamic dependencies were learned sequentially as data from 
additional informative flight test maneuvers were used in the global nonlinear modeling. 
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Figure 32. Prediction of fuzzy deceleration maneuver (Card 5) aerodynamics with model developed from 
data obtained during Cards 3 and 4 using the OFL modeling technique. 



Figure 33. Prediction of fuzzy deceleration maneuver (Card 5) aerodynamics with model developed from 
data obtained during Cards 3 and 4 using the MOF modeling technique. 
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Figure 34. Prediction of Card 6 data using model from previous maneuvers using the OFL technique. 



Figure 35. Prediction of Card 6 data using model from previous maneuvers using the MOF technique. 

This process continued in a similar way with the addition of data from a wind-up/wind-down turn (Card 7). 
For the two spin entries (Cards 8 and 9), the updated models predicted the response well, except for the regions that 
clearly were not covered in the previous flight test maneuvers. The spin entries included angles of attack, sideslip 
angles, and body-axis rotation rates that were far outside those experienced in the maneuvers to that point. Figure 
36 shows flight data for the right spin entry maneuver. Card 8. Compared to flight data from a fuzzy deceleration 
maneuver. Fig. 6 for instance, it is apparent that the ranges of explanatory variables achieved during the spin entry 
maneuver were not encompassed within flight data from previous maneuvers. 
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Figure 36. Right spin entry, Card 8. 


The model predictions of the right spin entry maneuver are shown in Figs. 37 and 38 for OFL and MOF, 
respectively. Inspection of the time-history matches of the model and flight data show that before the conditions 
became extreme, the global models from both modeling techniques exhibited excellent predictions. Throughout the 
entire maneuver, the OFL model does well predicting the longitudinal coefficients through extrapolation. During 
the extreme parts of the maneuver, both methods predict poorly, but similar to the flight data trends. This is an 
example of the extrapolation capabilities of the two modeling methods to conditions outside the modeled 
explanatory variable subspace. Similar prediction results were obtained for the left spin entry maneuver. Card 9. 

With no further model updates, another SPAZ maneuver was flown. Card 13, in the opposite turning direction 
from Card 4. The predictions for that maneuver with the OFL and MOF models are shown in Figs. 39 and 40, 
respectively. The OFL method produced excellent predictions for all of the flight data, and the MOF model 
produced slightly degraded predictions. 
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Figure 37. Prediction of right spin entry maneuver (Card 8) using model from previous maneuvers 

using the OFL technique. 



Figure 38. Prediction of right spin entry maneuver (Card 8) using model from previous maneuvers 

using the MOF technique. 
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Figure 39. Prediction of SPAZ maneuver (Card 13) using model from previous maneuvers using the OFL 

technique. 



Figure 40. Prediction of SPAZ maneuver (Card 13) using model from previous maneuvers using the MOF 

technique. 

Next, the models developed up through Card 13 were used to predict responses for some selected flight 
conditions in the flight envelope. Short, large amplitude, fuzzy maneuvers were flown at various flight speeds from 
cruise to near stall, and the predictions of the aerodynamics using the identified global models were compared to 
aerodynamics derived from flight data (Figs. 41 through 50). The predictions were generally very good with both 
modeling techniques, down to 100 KIAS. Although the model prediction for typically looks reasonable, the 
modeling often fails to meet the desired model metric targets. As discussed previously, this is likely due to 
problems in dealing with the thrust effects and A-axis aerodynamics in a combined manner, and may be indicative of 
missing needed explanatory variables, interaction effects, or issues with accuracy of measurements being used to 
generate the models. 
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Figures 41 and 42 show the model predictions for a fiizzy maneuver at 250 KIAS. The OFL and MOF models 
exhibit excellent predictions, except for the model, as noted earlier. This could be indicative of a need for more 
data in the identification process or possibly more degrees of freedom in the modeling, which for OFL would be 
more cells, and for MOF would be more spline functions available to MOF in the pool of candidate modeling terms. 



Figure 41. Prediction of maneuver at 250 KIAS using the OFL model. 



Figure 42. Prediction of maneuver at 250 KIAS using the MOF model. 
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Similar model prediction results were seen for maneuvers at 200 KIAS, as shown in Figs. 43 and 44. 




Figure 44. Prediction of maneuver at 200 KIAS using the MOF model. 


49 

American Institute of Aeronautics and Astronautics 












The model predictions for 150 KIAS are shown in Figs. 45 and 46. Prediction for both modeling methods is again 
very good. 



Figure 45. Prediction of maneuver at 150 KIAS using the OFL model. 



Figure 46. Prediction of maneuver at 150 KIAS using the MOF model. 


At 100 KIAS, Figs. 47 and 48 show degraded prediction for normal force coefficient . This may be an 
indicator of more data needed, or issues with instrumentation accuracy or repeatability, but in general, both 
modeling methods gave similar model predictions. A set of pitch, yaw, and roll doublets were also performed at 
100 KIAS. The prediction results for both modeling methods are shown in Figs. 49 and 50. The predictions are 
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very good, with the OFL model predicting all components except within the modeling metric targets. The MOF 
model prediction did very well with the lateral/directional coefficients, but showed more error in and . 



Figure 47. Prediction of maneuver at 100 KIAS using the OFL model. 



Figure 48. Prediction of maneuver at 100 KIAS using the MOF model. 


51 

American Institute of Aeronautics and Astronautics 












Figure 49. Prediction of doublets at 100 KIAS using the OFL model. 








tasi 






Figure 50. Prediction of doublets at 100 KIAS using the MOF model. 
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VI. Observations and Lessons Learned 

Several novel maneuvers were invented during this work, and all of them used the concept of simultaneous, 
multi-axis fuzzy excitation implemented by a pilot. Although the majority of the flight test maneuvers were flown 
by the lead author of this paper, two other test pilots from the NTPS were briefed on the maneuvers, and were 
provided with in-flight demonstrations of the new flight test maneuvers. All pilots flew the maneuvers very well. 
For example, the maneuvers shown in Figs. 19 through 31 were flown by a NTPS test pilot. This demonstrated that 
these new, highly efficient, global nonlinear modeling maneuvers are within the capabilities of a good test pilot. 

When working with global, nonlinear, high-dimensional models, the associated explanatory variable subspace is 
very complex. It is very difficult to determine when the explanatory variable subspace has been adequately covered 
by flight test maneuvers, and because the effects of any explanatory variable can be dependent on the values of other 
explanatory variables, the requirements for coverage are extensive. The approach taken in this work was to develop 
maneuvers that exercise a wide range of the explanatory variables as independently as possible, then test with 
different maneuvers to verify that the model predictions are adequate. Where deficiencies existed, more data were 
supplied to the modeling algorithms and the global nonlinear models were updated. Of course, having flight data 
and model predictions available in flight allows for real-time knowledge of the model adequacy, and indicates where 
more flight data are needed. This was a significant advantage for the general approach taken in this work. 

During flight 11, some basic aerobatic maneuvers were flown, including an inside loop, Cuban eight, 
Immelmann, and a clover leaf These are very smooth and slowly- varying maneuvers. For these maneuvers, the 
aerodynamic model identification for the moment coefficients was not as good as expected. Because the moment 
coefficient magnitudes were very small compared to the identification maneuvers used previously, it is possible that 
the lower values (typically by an order of magnitude) resulted in data with signal-to-noise ratios not well suited to 
modeling. To fully understand this issue, a more detailed characterization of onboard instrumentation system would 
be required. Some known problems that could be relevant are: friction in the angle of attack and angle of sideslip 
flow angle vanes, and variations in the air data with sideslip angle, which was handled in the onboard data 
conditioning by applying global Fourier smoothing to the recorded data. 

A comparison of the MOF and OFL modeling approaches was made by using both methods concurrently for 
each maneuver throughout the flight testing. In general, both approaches produced excellent global nonlinear 
models very rapidly, and were similar in performance. Note that although all of the modeling was done onboard the 
aircraft, the real-time maneuver monitoring and flight data recording was done at separate times from the global 
nonlinear modeling and prediction testing, because of limitations in programming the onboard tablet computer. This 
did not impact efficiency significantly, but could be improved upon by taking advantage of parallel processing, for 
example. 

In general, both modeling methods had the most difficulty with Fx modeling. The difficulty with Fx 
modeling was likely the combined effect of a strongly nonlinear aerodynamic component Cx that is generally 
difficult to model from flight data, combined with jet engine thrust needing a very different model form, and 
significant interaction effects between these two components. All of that was being characterized with a single 
global nonlinear Fx model expected to be valid for large portions of the flight envelope. 

As would be expected, each of the two methods used for onboard global nonlinear aerodynamic modeling has 
advantages and disadvantages. The fuzzy logic approach is a very powerful nonlinear modeling technique, but to 
implement a completely unrestrained fuzzy logic modeling approach, the current formulation requires much more 
computational time to develop the number of membership functions and fuzzy cells needed to adequately model 
each maneuver. In the current application, a candidate number of membership functions for each explanatory 
variable was pre-defmed based on previous experience. Also, in the fully unconstrained fuzzy modeling 
formulation, the addition of data to update a model requires keeping the old data instead of allowing for a recursive 
approach. This limitation of having to maintain the entire history of data might be mitigated with further 
development. The approach taken with the two modeling techniques is different as well. The MOF approach 
identifies a global model directly, using polynomials and splines with orthogonality concepts to identify a single 
model valid over the full range of the explanatory variables. The OFL approach splits the explanatory variable data 
into sub-ranges and ties together many local models using fuzzy logic, to create a global nonlinear model. For that 
reason, the OFL modeling can implement local model updates without degrading the rest of the model, whereas 
MOF does not have that capability in its current form. However MOF models typically use many fewer model 
parameters than OFL models, and therefore have easier physical interpretation, and typically extrapolate better 
outside the range of modeling data. 
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VII. Conclusions 

Research was conducted to develop efficient flight test techniques and onboard nonlinear modeling methods for 
identifying global aerodynamic models for aircraft from flight data. Flight data were obtained using piloted multi- 
axis perturbation inputs applied during slowly-varying flight conditions, resulting in rich data for global nonlinear 
aerodynamic modeling. This novel maneuver type exhibits a combination of low pair-wise correlations among 
explanatory variables, multi-axis excitation, and coverage of a large portion of the explanatory variable subspace for 
aerodynamic modeling, corresponding to large portions of the aircraft flight envelope. These characteristics make 
the maneuver very effective and efficient for global aerodynamic modeling. For example, using just one multi-axis 
fuzzy powered deceleration maneuver, a global nonlinear aerodynamic model valid for most of the normal operating 
range of the aircraft was identified from flight data alone. 

Two nonlinear modeling techniques were applied to identify global aerodynamic models - one using fuzzy 
logic, and the other using multivariate orthogonal functions. Both methods achieved excellent model fits to the data 
and excellent model prediction results, all in real time onboard the aircraft. The resulting global aerodynamic 
models were incorporated into a nonlinear simulation that exhibited good prediction capability for data from flight 
maneuvers that were not used in the modeling process. 

Flight test results show promise of substantial improvements in the capability to identify global aerodynamic 
models in flight regimes that traditional methods cannot even obtain, and in the rapid development of nonlinear 
models over large flight envelopes using very little flight test time. Maturation of this technology could result in 
substantial cost and time savings in flight test programs, and enable self-learning aircraft, resulting in reduced 
development cost for new aircraft, and more robust and safe flight operations. 

Piloted global modeling flight test maneuvers were used to collect data over relatively large ranges of the 
explanatory variables for all six rigid-body degrees of freedom of the aircraft simultaneously. This resulted in 
efficient collection of highly informative dynamic flight data with decorrelated explanatory variables, for good 
model identifiability and modeling accuracy. The nonlinear aerodynamic modeling methods identified accurate 
models without any specification of the model structures, resulting in a real-time onboard global aerodynamic 
modeling capability. This capability was developed and refined so that it was possible to identify global models for 
all six rigid-body degrees of freedom, including engine thrust modeling, for three aircraft configurations (cruise: 
flaps up, gear up, take-off: flaps half, gear up, and approach/landing: flaps full, gear down), using two different 
global modeling techniques, and including validation of all identified models, all in a single test flight lasting less 
than 1 hour. 

Accurate global models identified from flight data alone could be used in rapid and relatively low-budget 
unmanned aircraft programs to save development time and money by making it possible to generate an accurate 
nonlinear aircraft simulation from a single flight, without the need for extensive wind tunnel testing or aerodynamic 
calculations. The capability also has important implications for aircraft safety, because the technique could be used 
to generate an onboard model of the global aerodynamics for an aircraft. The approach could be used to account for 
the particular geometry or flight environment of an individual aircraft, as well as provide a capability to monitor and 
account for changes in the aircraft due to failures, damage, and airframe icing, for example. A global aerodynamic 
model identified from flight data alone could also be compared to aerodynamic databases generated from wind 
tunnel data and aerodynamic calculations to improve the fidelity of the ground-based aerodynamic prediction 
methods. This could be done through identifying where there are significant differences between the database 
generated using ground-based methods and the global model identified from flight data alone, as well as where the 
models agreed well. Such comparisons could also be used to estimate full-scale Reynolds number effects, or 
artifacts of wind tunnel testing, such as sting interference and wall effects. 

The flight testing and global modeling techniques developed in this work could also be applied to other aircraft, 
such as regional jets and turboprop aircraft. The improved efficiency makes the idea of identifying flight 
simulations from flight data for each particular aircraft (not just every aircraft of a particular type and model) a 
reasonable proposition. The process also has great potential for time and cost savings in certifying flight 
simulations, because the flight test match criteria would be satisfied automatically by the fact that the global 
aerodynamic model was identified directly from flight data. 
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